General project set-up

# Get all libraries and sources required to run the script
        library(plyr)
        library(dplyr)
        library(reshape2)

        library(ggplot2)
        library(ggthemes)
        library(scales)

        library(lme4)
        library(lmerTest)
        library(emmeans)
    
        library(skimr)

# Default ggplot settings

    Fill.colour<-scale_colour_manual(values = c("black", "gray70", "gray35"))

    ggthe_bw<-theme(plot.background=element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          legend.box.background = element_rect(),
          panel.background =element_rect(fill = NA, color = "black")
          )+
    theme_bw()

1. DATA Exploration

1. Select qPCR info and define factors

# Import the data
  qPCR.variables  <- read.csv("Outputs/Ssid_SH_cell_ratio.csv", header = T)

# Remove blastate samples
  qPCR.blastate<-qPCR.variables[((qPCR.variables$Date=="2018-02-05") |
                                   (qPCR.variables$Date=="2018-03-08")), ]
  qPCR.variables <- droplevels(qPCR.variables[!rownames(qPCR.variables) %in% rownames(qPCR.blastate), ])  
  
 
# Variable types 
str(qPCR.variables)
## 'data.frame':    644 obs. of  17 variables:
##  $ Treatment : Factor w/ 3 levels "A","N","N+P": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Replicate : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Fragment  : Factor w/ 162 levels "Ss_20-1","Ss_20-11",..: 27 30 31 48 49 52 55 58 97 100 ...
##  $ Genotype  : Factor w/ 7 levels "Ss_20","Ss_22",..: 2 2 2 2 3 3 3 3 5 5 ...
##  $ Date      : Factor w/ 5 levels "2017-11-15","2017-12-14",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Days      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Time_Point: Factor w/ 5 levels "T13","T19","T22",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Phase     : Factor w/ 3 levels "Baseline","Heat",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Sample    : Factor w/ 644 levels "Ss_20-1_T5","Ss_20-1_T9",..: 109 124 128 188 191 204 217 232 380 393 ...
##  $ C.Ssid    : num  0 0 0 0 0 ...
##  $ D.Ssid    : num  0.0266 0.034 0.0239 0.0281 1.0265 ...
##  $ TotalSH   : num  0.0266 0.034 0.0239 0.0281 1.0265 ...
##  $ logC.SH   : num  NA NA NA NA NA ...
##  $ logD.SH   : num  -1.5756 -1.4689 -1.621 -1.5515 0.0114 ...
##  $ logSH     : num  -1.5756 -1.4689 -1.621 -1.5515 0.0114 ...
##  $ D.Prp     : num  1 1 1 1 1 ...
##  $ Community : Factor w/ 3 levels "C1","C3","D": 3 3 3 3 3 3 3 3 3 3 ...
  qPCR.variables$Genotype<-factor(qPCR.variables$Genotype, 
                                   levels=c("Ss_22","Ss_23","Ss_27", "Ss_28",
                                          "Ss_20", "Ss_24","Ss_30"))
  qPCR.variables$DaysF<-as.factor(qPCR.variables$Days)
  
  summary(qPCR.variables)
##  Treatment Replicate     Fragment    Genotype          Date    
##  A  :220   R1:333    Ss_20-17:  5   Ss_22:92   2017-11-15:140  
##  N  :228   R2:311    Ss_20-19:  5   Ss_23:96   2017-12-14:154  
##  N+P:196             Ss_20-22:  5   Ss_27:96   2018-02-01:151  
##                      Ss_20-26:  5   Ss_28:82   2018-02-23:103  
##                      Ss_20-29:  5   Ss_20:97   2018-03-06: 96  
##                      Ss_20-32:  5   Ss_24:85                   
##                      (Other) :614   Ss_30:96                   
##       Days        Time_Point       Phase              Sample   
##  Min.   :  0.00   T13:151    Baseline :140   Ss_20-1_T5  :  1  
##  1st Qu.: 29.00   T19:103    Heat     :199   Ss_20-1_T9  :  1  
##  Median : 78.00   T22: 96    Nutrients:305   Ss_20-11_T13:  1  
##  Mean   : 57.76   T5 :140                    Ss_20-11_T19:  1  
##  3rd Qu.:100.00   T9 :154                    Ss_20-11_T22:  1  
##  Max.   :111.00                              Ss_20-11_T9 :  1  
##                                              (Other)     :638  
##      C.Ssid             D.Ssid            TotalSH         
##  Min.   :0.000000   Min.   :0.000000   Min.   :0.0000034  
##  1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.0093391  
##  Median :0.002373   Median :0.006613   Median :0.0313298  
##  Mean   :0.046419   Mean   :0.059660   Mean   :0.1060796  
##  3rd Qu.:0.027225   3rd Qu.:0.044028   3rd Qu.:0.1286860  
##  Max.   :1.155796   Max.   :1.394495   Max.   :1.3944953  
##                                                           
##     logC.SH            logD.SH            logSH             D.Prp       
##  Min.   :-5.47264   Min.   :-4.7761   Min.   :-5.4726   Min.   :0.0000  
##  1st Qu.:-2.80684   1st Qu.:-2.2410   1st Qu.:-2.0297   1st Qu.:0.0000  
##  Median :-2.09749   Median :-1.6835   Median :-1.5040   Median :0.6792  
##  Mean   :-2.08904   Mean   :-1.7359   Mean   :-1.5221   Mean   :0.5242  
##  3rd Qu.:-1.27217   3rd Qu.:-1.0631   3rd Qu.:-0.8905   3rd Qu.:1.0000  
##  Max.   : 0.06288   Max.   : 0.1444   Max.   : 0.1444   Max.   :1.0000  
##  NA's   :184        NA's   :205                                         
##  Community DaysF    
##  C1:132    0  :140  
##  C3:171    29 :154  
##  D :341    78 :151  
##            100:103  
##            111: 96  
##                     
## 
  skim(qPCR.variables)
Data summary
Name qPCR.variables
Number of rows 644
Number of columns 18
_______________________
Column type frequency:
factor 10
numeric 8
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Treatment 0 1 FALSE 3 N: 228, A: 220, N+P: 196
Replicate 0 1 FALSE 2 R1: 333, R2: 311
Fragment 0 1 FALSE 162 Ss_: 5, Ss_: 5, Ss_: 5, Ss_: 5
Genotype 0 1 FALSE 7 Ss_: 97, Ss_: 96, Ss_: 96, Ss_: 96
Date 0 1 FALSE 5 201: 154, 201: 151, 201: 140, 201: 103
Time_Point 0 1 FALSE 5 T9: 154, T13: 151, T5: 140, T19: 103
Phase 0 1 FALSE 3 Nut: 305, Hea: 199, Bas: 140
Sample 0 1 FALSE 644 Ss_: 1, Ss_: 1, Ss_: 1, Ss_: 1
Community 0 1 FALSE 3 D: 341, C3: 171, C1: 132
DaysF 0 1 FALSE 5 29: 154, 78: 151, 0: 140, 100: 103

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Days 0 1.00 57.76 41.59 0.00 29.00 78.00 100.00 111.00 ▆▆▁▆▇
C.Ssid 0 1.00 0.05 0.12 0.00 0.00 0.00 0.03 1.16 ▇▁▁▁▁
D.Ssid 0 1.00 0.06 0.14 0.00 0.00 0.01 0.04 1.39 ▇▁▁▁▁
TotalSH 0 1.00 0.11 0.17 0.00 0.01 0.03 0.13 1.39 ▇▁▁▁▁
logC.SH 184 0.71 -2.09 1.08 -5.47 -2.81 -2.10 -1.27 0.06 ▁▃▇▇▅
logD.SH 205 0.68 -1.74 0.94 -4.78 -2.24 -1.68 -1.06 0.14 ▁▂▆▇▃
logSH 0 1.00 -1.52 0.83 -5.47 -2.03 -1.50 -0.89 0.14 ▁▁▃▇▅
D.Prp 0 1.00 0.52 0.47 0.00 0.00 0.68 1.00 1.00 ▇▁▁▁▇

2. Exploratory graphs

  • Histograms
  HistoL_SH<-qplot(logSH, data=qPCR.variables, binwidth=0.15)
   HistoL_SH + facet_grid(Treatment~Date)

  ggplot(qPCR.variables, aes(logSH, fill = Treatment , colour = Treatment)) +
   geom_density(alpha = 0.1) + facet_wrap(~Date) + ggthe_bw

  • Log SH
logSHGenotype <- ggplot(qPCR.variables, aes(Days, logSH)) +
  geom_line(aes(colour=Fragment))+geom_point(aes(shape=factor(Community), colour=factor(Community)))+
  # geom_jitter(aes(colour=factor(Replicate))) +
      # stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
      # stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
      facet_grid(Treatment~Genotype) +
      ggthe_bw +theme(legend.position = "none" )

logSHGenotype + ylab("Relative log10 (S:H)") + xlab("Treatment") +  
      theme(axis.title.y=element_text(size=12), legend.position="none")

logSHTreatment <- ggplot(qPCR.variables, aes(Date, logSH, colour=factor(Treatment))) +
       # geom_jitter(aes(colour=factor(Treatment))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
       stat_summary(fun.y=mean, geom="line", size =1) +
       theme_bw() 
logSHTreatment + facet_wrap(~Genotype)

logSH_Replicate<- ggplot(qPCR.variables, aes (Days, logSH, 
                                             colour=factor(Replicate))) +
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
  stat_summary(fun.y=mean, geom="line") + facet_grid (~Treatment) + ggthe_bw
logSH_Replicate 

logSH_Replicate+ facet_grid(~Community)

logSHTreatment <- ggplot (qPCR.variables, aes(Days, logSH, colour=Treatment)) + 
  ggtitle("A.") +  # geom_point(alpha=0.3) +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 1)) +
  stat_summary(fun.y=mean, geom="line", position = position_dodge(width = 5)) + 

  scale_y_continuous(breaks = seq(-5, 0.5, 0.5),
                    expand = c(0,0),
                  name=("log 10 (Total S/H cell ratio")) +
  scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0)) +  
 
      annotate("segment", x = 2, xend = 91, y = -3.5, yend = -3.5,
                colour = "gray90")+
      annotate("segment", x = 79, xend = 90, y = -3.5, yend = -3,
                colour = "gray90")+
      annotate("segment", x = 91, xend = 110, y = -3, yend = -3,
                colour = "gray90")
logSHTreatment

#logSHTreatment + facet_wrap(~Genotype)
logSHTreatment + facet_wrap(~Community)

  • Log C
logC_Treatment <- ggplot(qPCR.variables,
                         aes(Days, logC.SH, colour=factor(Fragment))) +
     geom_line()+ 
     #geom_jitter(aes(colour=factor(Replicate))) +
         #stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
         #stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
         ggthe_bw +  theme(legend.position = "none" )
logC_Treatment

logC_Treatment + facet_grid(Genotype~Treatment)

logC_Treatment + facet_grid(Community~Treatment)

  • Log D
logD_genotype <- ggplot(qPCR.variables,
                         aes(Days, logD.SH, colour=factor(Fragment))) +
     geom_line()+ 
     #geom_jitter(aes(colour=factor(Replicate))) +
         #stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
         #stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
         ggthe_bw +  theme(legend.position = "none" )
logD_genotype

logD_genotype + facet_grid(Genotype~Treatment)

logD_Treat <- ggplot(qPCR.variables,
                         aes(Days, logD.SH, colour=factor(Treatment))) +
     #geom_jitter(aes(colour=factor(Replicate))) +
        stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
        stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
        ggthe_bw +  theme(legend.position = "none" )
logD_Treat

logD_Treat + facet_grid(Community~Treatment)

*D proportion

D.PTreatment <- ggplot(qPCR.variables, aes(Days, D.Prp, colour=Fragment)) +
        ggthe_bw +  theme(legend.position="none") +  
        geom_jitter(aes(shape=factor(Replicate)), alpha=0.5) +
        geom_line() + facet_grid (Genotype~Treatment) +
  
      annotate("segment", x = 2, xend = 91, y = 0, yend = 0,
                colour = "gray90")+
      #annotate("text", x = 45, y = 0.2, label = "Nutrients", size=3)+
    
      annotate("segment", x = 79, xend = 90, y = 0, yend = 0.3,
                colour = "gray90")+
    
      annotate("segment", x = 91, xend = 110, y = 0.3, yend = 0.3,
                colour = "gray90")+
      #annotate("text", x = 99, y = 0.3, label = "Heat", size=3)+

      scale_y_continuous(breaks = seq(0, 1, 0.3),
                     name=("Durusdinium proportion (D/H)/(S/H)")) +
      scale_x_continuous(name="Days in the experiment",
                           limits = c(-2,113),
                         breaks = seq(0, 110, 45),  
                         expand = c(0, 0))
D.PTreatment

2. Subset data

# Removing weird points and colonies not used

  qPCR.variables_2<-subset(qPCR.variables, Date!="2017-12-14")
  qPCR.variables_2<-subset(qPCR.variables_2, TotalSH<0.8)
  
  SH.0<-subset(qPCR.variables_2, Days<2)
  SH.C<-subset(qPCR.variables_2, Days<79)
    SH.C<-subset(SH.C, Days>77)
  SH.H<-subset(qPCR.variables_2, Days>80)
 
  
logSH_Genotype<- ggplot(qPCR.variables_2, aes (Days, logSH, 
                     colour=factor(Fragment))) +
  theme_bw() + geom_line()+ facet_grid (Genotype~Treatment)#+
logSH_Genotype + theme(legend.position = "none" ) + geom_point()

Log S/H graph

logSHTreatment <- ggplot (qPCR.variables_2, aes(Days, logSH, colour=Treatment)) +  theme_bw() + ggtitle("A.")+        Fill.colour+ ggthe_bw+
    theme(plot.background=element_blank(), 
            panel.border = element_blank(),
            axis.line = element_line(colour = "black"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position="bottom",
            strip.background = element_rect(fill="white")) +
       #geom_jitter(aes(colour=factor(Replicate))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 5)) +
       stat_summary(fun.y=mean, geom="line") + 
      
  annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +

      scale_y_continuous(breaks = seq(-5, 0.3, 0.5),
                    expand = c(0,0),
                   name=("Symbiodiniaceae / S. siderea (S/H) cell ratio")) +
      scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0))
  logSHTreatment 

  logSHTreatment + facet_grid(~Community)

  #ggsave(file="5.3A_AH.svg", plot=logSHTreatment, width=3.0, height=6)
  #ggsave(file="5.3A_AH_NoGenotype.svg", plot=logSHTreatment, width=3.0, height=3)

Log C/H graph

CHTreatment_2 <- ggplot (qPCR.variables_2, aes(Days, logC.SH, colour=Treatment)) +  theme_bw() + ggtitle("C.")+ Fill.colour+ ggthe_bw+
  annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
       
       #geom_jitter(aes(colour=factor(Replicate))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 5)) +
       stat_summary(fun.y=mean, geom="line") +
      scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0))
  CHTreatment_2  

  #CHTreatment_2 + facet_grid (Genotype~Community)

  CHTreatment_2 + facet_grid(~Community)

  #ggsave(file="5.3A_AH_2.svg", plot=SHTreatment_2, width=3.0, height=6)
  #ggsave(file="5.3C_SH_NoGenotype_2.svg", plot=SHTreatment_2, width=3.0, height=3)

Log D/H graph

DHTreatment_2 <- ggplot (qPCR.variables_2, aes(Days, logD.SH, colour=Treatment)) +  theme_bw() + ggtitle("C.")+ Fill.colour+ ggthe_bw+
  annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
       
       #geom_jitter(aes(colour=factor(Replicate))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 5)) +
       stat_summary(fun.y=mean, geom="line") +
      scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0))
  DHTreatment_2 + facet_grid (Genotype~Community)

  DHTreatment_2 + facet_grid(~Community)

D Proportion

D.PTreatment <- ggplot(qPCR.variables, aes(Days, D.Prp, colour=Fragment)) +
        ggthe_bw +  theme(legend.position="none") +  
        geom_jitter(aes(shape=factor(Treatment)), alpha=0.5) +
        geom_line()+ 
      annotate("segment", x = 2, xend = 91, y = 0, yend = 0,
                colour = "gray90")+
      #annotate("text", x = 45, y = 0.2, label = "Nutrients", size=3)+
    
      annotate("segment", x = 79, xend = 90, y = 0, yend = 0.3,
                colour = "gray90")+
    
      annotate("segment", x = 91, xend = 110, y = 0.3, yend = 0.3,
                colour = "gray90")+
      #annotate("text", x = 99, y = 0.3, label = "Heat", size=3)+

      scale_y_continuous(breaks = seq(0, 1, 0.3),
                     name=("Durusdinium proportion (D/H)/(S/H)")) +
      scale_x_continuous(name="Days in the experiment",
                           limits = c(-2,113),
                         breaks = seq(0, 110, 45),  
                         expand = c(0, 0))
      D.PTreatment + facet_grid (Genotype~Replicate)

D.PTreatment

3. GLM models

Baseline

LME_BaseLine<-lmer(logSH ~ Treatment + Community + (1|Genotype) + (1|Replicate), 
                     data=SH.0)
   step(LME_BaseLine)
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        8  -75.016 166.03                         
## (1 | Genotype)           0    7  -82.009 178.02 13.986  1  0.0001842 ***
## (1 | Replicate)          0    7 -109.286 232.57 68.539  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated  Sum Sq Mean Sq NumDF   DenDF F value Pr(>F)
## Treatment          1 0.59627 0.29813     2 126.566  2.0641 0.1312
## Community          2 0.79889 0.39945     2   6.535  2.7209 0.1381
## 
## Model found:
## logSH ~ (1 | Genotype) + (1 | Replicate)
   anova(LME_BaseLine) # Treatments is not significant
   ranova(LME_BaseLine) # Treatments is significant!
   LME_BaseLine<-lmer(logSH ~ Genotype + (1|Replicate), 
                     data=SH.0)
   step(LME_BaseLine)
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        9  -72.767 163.53                         
## (1 | Replicate)          0    8 -106.895 229.79 68.256  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##          Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Genotype          0 8.6374  1.4396     6   129  9.7555 7.505e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Genotype + (1 | Replicate)
   anova(LME_BaseLine) # Treatments is not significant
   ranova(LME_BaseLine) # Treatments is significant!
# Multicomp
Ssid.SH.C.emm<-emmeans(LME_BaseLine, ~Genotype)
      #contrast(Ofav.SH.C.emm, "tukey")
    Ssid.SH.C_groups<-cld(Ssid.SH.C.emm, by=NULL) # compact-letter display
    Ssid.SH.C_groups
# Effect plot
  plot(emmeans(LME_BaseLine, ~Genotype), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()

logSHTreatmentBL <- ggplot (SH.0, aes(Genotype, logSH, colour=Genotype)) + 
  #ggtitle("A.") + 
  ggthe_bw +  theme(legend.position="bottom") +
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentBL

logSHTreatmentBL + facet_grid(~Community)

C: Effect of nutrient treatments at control temperature

LME_Ssid.C<-lmer(logSH ~ Treatment + Community + (1|Replicate) + (1|Genotype), 
                     data=SH.C)
     step (LME_Ssid.C) #  Remove Genotype
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        8  -73.147 162.29                         
## (1 | Genotype)           1    7  -73.264 160.53  0.235  1     0.6278    
## (1 | Replicate)          0    6 -105.884 223.77 65.240  1   6.63e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment          1 0.7873  0.3936     2   144  2.8627   0.06037 .  
## Community          0 9.4779  4.7390     2   146 33.6053 9.879e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Community + (1 | Replicate)
     anova (LME_Ssid.C)
     ranova (LME_Ssid.C)
LME_Ssid.C<-lmer(logSH ~ Treatment + Community + (1|Replicate),
                     data=SH.C)
     step (LME_Ssid.C) #  Remove replicate
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC   LRT Df Pr(>Chisq)    
## <none>                        7  -73.264 160.53                        
## (1 | Replicate)          0    6 -105.884 223.77 65.24  1   6.63e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment          1 0.7873  0.3936     2   144  2.8627   0.06037 .  
## Community          0 9.4779  4.7390     2   146 33.6053 9.879e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Community + (1 | Replicate)
     anova (LME_Ssid.C)
     ranova (LME_Ssid.C)
# Multicomp
Ssid.SH.C.emm<-emmeans(LME_Ssid.C, ~Community)
      #contrast(Ssid.SH.C.emm, "tukey")
    Ssid.SH.C_groups<-cld(Ssid.SH.C.emm, by=NULL) # compact-letter display
    Ssid.SH.C_groups
# Effect plot
  plot(emmeans(LME_Ssid.C, ~Treatment), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()

logSHTreatmentC <- ggplot (SH.C, aes(Treatment, logSH, colour=Treatment, shape=Treatment)) + 
  ggtitle("A.")+ Fill.colour +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentC

logSHGenotypeC <- ggplot (SH.C, aes(Genotype, logSH, colour=Genotype)) + 
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHGenotypeC

H: Effect of pre-exposure to nutrient treatments during heat challenge

LME_Ssid.H<-lmer(logSH ~ Treatment + Community + DaysF+ (1|Replicate) + (1|Genotype), 
                     data=SH.H)
     step (LME_Ssid.H) #  Remove Replicate
## Backward reduced random-effect table:
## 
##                 Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        9 -150.52 319.04                         
## (1 | Replicate)          1    8 -150.52 317.04  0.000  1     0.9999    
## (1 | Genotype)           0    7 -178.47 370.93 55.896  1  7.641e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Treatment          0 2.8371  1.4186     2 186.93  6.4838  0.001895 ** 
## Community          0 5.8196  2.9098     2 189.17 13.2996 3.942e-06 ***
## DaysF              0 9.3166  9.3166     1 187.32 42.5827 6.161e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Community + DaysF + (1 | Genotype)
     anova (LME_Ssid.H)
     ranova (LME_Ssid.H)
LME_Ssid.H<-lmer(logSH ~ Treatment + Community + DaysF+ (1|Genotype),
                     data=SH.H)
     step (LME_Ssid.H) #  Remove replicate
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                       8 -150.52 317.04                         
## (1 | Genotype)          0    7 -178.47 370.93 55.896  1  7.641e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Treatment          0 2.8371  1.4186     2 186.93  6.4838  0.001895 ** 
## Community          0 5.8196  2.9098     2 189.17 13.2996 3.942e-06 ***
## DaysF              0 9.3166  9.3166     1 187.32 42.5827 6.161e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Community + DaysF + (1 | Genotype)
     anova (LME_Ssid.H)
     ranova (LME_Ssid.H)
# Multicomp
Ssid.SH.H.emm<-emmeans(LME_Ssid.H, ~Treatment)
      #contrast(Ssid.SH.H.emm, "tukey")
    Ssid.SH.H_groups<-cld(Ssid.SH.H.emm, by=NULL) # compact-letter display
    Ssid.SH.H_groups
# Effect plot
  plot(emmeans(LME_Ssid.H, ~Treatment), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()

logSHTreatmentH <- ggplot (SH.H, aes(Treatment, logSH, colour=Treatment, shape=Treatment)) + 
  ggtitle("A.")+ Fill.colour +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-3, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentH

logSHTreatmentH + facet_grid(~Community)

logSHTreatmentH + facet_grid(DaysF~Community) #+  geom_jitter()

logSHGenotypeH <- ggplot (SH.H, aes(Community, logSH, colour=Genotype)) + 
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHGenotypeH

logSHGenotypeH + facet_grid(~DaysF)

All phases

  • Dasys as a factor
LM_Ssid.qpCR<-lmer(logSH ~ Treatment * DaysF * Community + (1|Replicate) + (1|Genotype/Fragment), 
                   data= qPCR.variables_2)
    step(LM_Ssid.qpCR)
## Backward reduced random-effect table:
## 
##                         Eliminated npar  logLik    AIC    LRT Df
## <none>                               40 -345.92 771.83          
## (1 | Replicate)                  1   39 -345.92 769.83  0.000  1
## (1 | Fragment:Genotype)          2   38 -345.92 767.83  0.000  1
## (1 | Genotype)                   0   37 -376.42 826.84 61.013  1
##                         Pr(>Chisq)    
## <none>                                
## (1 | Replicate)                  1    
## (1 | Fragment:Genotype)          1    
## (1 | Genotype)           5.669e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:DaysF:Community          1  2.3178  0.1932    12 443.35  0.9151
## Treatment:Community                2  0.7151  0.1788     4 455.41  0.8489
## Treatment:DaysF                    0  5.3896  0.8983     6 459.60  4.2710
## DaysF:Community                    0 25.2283  4.2047     6 459.70 19.9919
##                              Pr(>F)    
## Treatment:DaysF:Community  0.531574    
## Treatment:Community        0.494734    
## Treatment:DaysF            0.000335 ***
## DaysF:Community           < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF + 
##     DaysF:Community
    anova(LM_Ssid.qpCR)
    ranova(LM_Ssid.qpCR)
LM_Ssid.qpCR<-lmer(logSH ~ Treatment + DaysF + Community + 
                     (1 | Genotype) + 
                     Treatment:DaysF +  DaysF:Community,
                   data= qPCR.variables_2)
      step(LM_Ssid.qpCR)
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      22 -344.90 733.79                         
## (1 | Genotype)          0   21 -376.15 794.30 62.504  1  2.659e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                 Eliminated  Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)
## Treatment:DaysF          0  5.3896  0.8983     6 459.6   4.271  0.000335
## DaysF:Community          0 25.2283  4.2047     6 459.7  19.992 < 2.2e-16
##                    
## Treatment:DaysF ***
## DaysF:Community ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF + 
##     DaysF:Community
      anova(LM_Ssid.qpCR)
      ranova(LM_Ssid.qpCR)
      summary(LM_Ssid.qpCR)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF +  
##     DaysF:Community
##    Data: qPCR.variables_2
## 
## REML criterion at convergence: 689.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3945 -0.6829 -0.0445  0.7223  3.3266 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Genotype (Intercept) 0.4366   0.6608  
##  Residual             0.2103   0.4586  
## Number of obs: 486, groups:  Genotype, 7
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            -1.91537    0.28599   8.61510  -6.697 0.000109 ***
## TreatmentN             -0.07283    0.09326 459.24915  -0.781 0.435249    
## TreatmentN+P           -0.23402    0.09857 459.28241  -2.374 0.018003 *  
## DaysF78                 0.18484    0.14532 459.31173   1.272 0.204047    
## DaysF100               -0.87819    0.15269 459.54363  -5.752 1.62e-08 ***
## DaysF111               -1.16397    0.15636 459.47375  -7.444 4.85e-13 ***
## CommunityC3             1.46403    0.22145 395.52840   6.611 1.24e-10 ***
## CommunityD              0.16239    0.14456 465.99482   1.123 0.261877    
## TreatmentN:DaysF78      0.09794    0.13091 459.24978   0.748 0.454730    
## TreatmentN+P:DaysF78    0.08522    0.13522 459.27412   0.630 0.528863    
## TreatmentN:DaysF100     0.48788    0.14188 459.32643   3.439 0.000638 ***
## TreatmentN+P:DaysF100   0.43816    0.15142 459.30234   2.894 0.003990 ** 
## TreatmentN:DaysF111     0.09205    0.14651 460.63193   0.628 0.530117    
## TreatmentN+P:DaysF111   0.48196    0.15444 459.53356   3.121 0.001917 ** 
## DaysF78:CommunityC3     0.35575    0.16130 459.31491   2.205 0.027917 *  
## DaysF100:CommunityC3   -0.12994    0.17289 459.43506  -0.752 0.452701    
## DaysF111:CommunityC3   -1.18260    0.19715 459.34228  -5.999 4.04e-09 ***
## DaysF78:CommunityD      0.10776    0.14580 459.32787   0.739 0.460215    
## DaysF100:CommunityD     0.02428    0.15623 459.39217   0.155 0.876552    
## DaysF111:CommunityD     0.31536    0.15824 460.63181   1.993 0.046859 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Emmeans
Ssid.YII.emm<-emmeans(LM_Ssid.qpCR, ~Treatment*DaysF*Community)
     #contrast(Ssid.YII.emm, "tukey")

# Pairwise comparisons
    #pairs(Ssid.YII.emm) # same than contrast(Ssid.YII.emm, "tukey")
  plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF), comparisons = TRUE) +
        theme_bw() # Tukey comparission (do not trust CI to compare EMMs)

    plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF|Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(Community~DaysF) +
      theme_bw()

  #CLD  
  SH_Multicomp<-cld(Ssid.YII.emm, by=NULL) # compact-letter display
      SH_Multicomp<-SH_Multicomp[order(SH_Multicomp$DaysF, 
                                       SH_Multicomp$Community,
                                       SH_Multicomp$Treatment), ]
    write.csv(SH_Multicomp, "Ssid_SH_Multicomp.csv")  
 # 2. Predict values:
    pred_Ssid1 <- predict(LM_Ssid.qpCR, re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1 <- bootMer(LM_Ssid.qpCR, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ssid1 - std.err*1.96
    CI.hi_1 <- pred_Ssid1 + std.err*1.96

  #Plot
  Model_Ss_1b_plot<- ggplot(
    qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    #scale_y_continuous(name=expression(~italic("Fv / Fm")),
    #                   limits = c(0.5, 0.61), 
    #                   breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    #scale_x_continuous("Days in the experiment", limits = c(0, 78),
    #                 breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                  linetype=1, alpha=1) + 
    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_1b_plot + facet_wrap(~Community)

  • Days as continous
LM_Ssid.qpCR2<-lmer(logSH ~ Treatment * Days * Community + (1|Replicate) + (1|Genotype/Fragment), 
                   data= qPCR.variables_2)
    step(LM_Ssid.qpCR2)
## Backward reduced random-effect table:
## 
##                         Eliminated npar  logLik    AIC    LRT Df
## <none>                               22 -567.31 1178.6          
## (1 | Replicate)                  1   21 -567.31 1176.6  0.000  1
## (1 | Fragment:Genotype)          2   20 -567.31 1174.6  0.000  1
## (1 | Genotype)                   0   19 -588.46 1214.9 42.303  1
##                         Pr(>Chisq)    
## <none>                                
## (1 | Replicate)                  1    
## (1 | Fragment:Genotype)          1    
## (1 | Genotype)           7.817e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                          Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:Days:Community          1  1.9970  0.4992     4 461.13  1.0674
## Treatment:Community               2  3.0175  0.7544     4 465.28  1.6124
## Days:Community                    3  2.0249  1.0125     2 469.82  2.1535
## Community                         0 23.5575 11.7787     2 371.44 24.9107
## Treatment:Days                    0  3.1263  1.5631     2 471.54  3.3058
##                             Pr(>F)    
## Treatment:Days:Community   0.37206    
## Treatment:Community        0.16997    
## Days:Community             0.11722    
## Community                7.047e-11 ***
## Treatment:Days             0.03752 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
    anova(LM_Ssid.qpCR2)
    ranova(LM_Ssid.qpCR2)
LM_Ssid.qpCR2<-lmer(logSH ~Treatment + Days + Community + (1|Genotype) + Treatment:Days,
                   data= qPCR.variables_2)
      step(LM_Ssid.qpCR2)
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      10 -542.25 1104.5                         
## (1 | Genotype)          0    9 -563.55 1145.1 42.618  1  6.656e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                Eliminated  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## Community               0 23.5575 11.7787     2 371.44 24.9107 7.047e-11
## Treatment:Days          0  3.1263  1.5631     2 471.54  3.3058   0.03752
##                   
## Community      ***
## Treatment:Days *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
      anova(LM_Ssid.qpCR2)
      ranova(LM_Ssid.qpCR2)
      summary(LM_Ssid.qpCR2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
##    Data: qPCR.variables_2
## 
## REML criterion at convergence: 1084.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5339 -0.6119 -0.1030  0.6531  2.3803 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Genotype (Intercept) 0.6894   0.8303  
##  Residual             0.4728   0.6876  
## Number of obs: 486, groups:  Genotype, 7
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -1.889571   0.359416   8.086576  -5.257 0.000741 ***
## TreatmentN         -0.097590   0.136063 471.087555  -0.717 0.473579    
## TreatmentN+P       -0.267215   0.143404 471.032500  -1.863 0.063032 .  
## Days               -0.008154   0.001207 471.192021  -6.755 4.21e-11 ***
## CommunityC3         1.865263   0.278743 295.864217   6.692 1.10e-10 ***
## CommunityD          0.305446   0.160114 461.670216   1.908 0.057053 .  
## TreatmentN:Days     0.003255   0.001699 471.935734   1.916 0.055934 .  
## TreatmentN+P:Days   0.004351   0.001797 471.117838   2.421 0.015863 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnN TrtN+P Days   CmmnC3 CmmntD TrtN:D
## TreatmentN  -0.183                                          
## TreatmntN+P -0.183  0.480                                   
## Days        -0.246  0.591  0.564                            
## CommunityC3 -0.356 -0.027  0.001  0.051                     
## CommunityD  -0.363 -0.009  0.006  0.044  0.563              
## TrtmntN:Dys  0.145 -0.834 -0.400 -0.705  0.070 -0.007       
## TrtmntN+P:D  0.137 -0.399 -0.838 -0.669  0.028  0.033  0.477
 # 2. Predict values:
    pred_Ssid1 <- predict(LM_Ssid.qpCR2,re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1 <- bootMer(LM_Ssid.qpCR2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ssid1 - std.err*1.96
    CI.hi_1 <- pred_Ssid1 + std.err*1.96

  #Plot
  Model_Ss_1b_plot<- ggplot(
    qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    #scale_y_continuous(name=expression(~italic("Fv / Fm")),
    #                   limits = c(0.5, 0.61), 
    #                   breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    #scale_x_continuous("Days in the experiment", limits = c(0, 78),
    #                 breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                  linetype=1, alpha=1) + 
    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_1b_plot + facet_wrap(Genotype~Community)